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Abstract 

Standard lattice fermion algorithms run into the well-known sign problem 
at real chemical potential. In this paper we investigate the possibility of 
using imaginary chemical potential, and argue that it has advantages over 
other methods, particularly for probing the physics at finite temperature as 
well as density. As a feasibility study, we present numerical results for the 
partition function of the two-dimensional Hubbard model with imaginary 
chemical potential. 

We also note that systems with a net imbalance of isospin may be 
simulated using a real chemical potential that couples to I3 without suffering 
from the sign problem. 



PACS numbers: 11.15.Ha, 71.10.Fd 



1 Introduction 



The behavior of fermions in the presence of a chemical potential is relevant to 
condensed matter physics (Hubbard model away from half-filling) and particle 
physics (high quark density systems such as the early universe, neutron stars, and 
heavy- ion collisions). Furthermore, a remarkably rich phase structure has been 
conjectured for QCD at finite temperature and density ||. 

The only reliable non-perturbative approach to QCD is the numerical Monte- 
Carlo evaluation of the functional integral using a lattice regulator. Unfortunately, 
standard Monte-Carlo methods become inapplicable at finite quark density, since 
in the presence of a real chemical potential the measure is no longer positive. One 
approach to this problem is the "Glasgow method" ||, in which the partition 
function is expanded in powers of e^, and the coefficients are evaluated by 
Monte-Carlo, using an ensemble of configurations weighted by the /i = action. 
Simulations using this method have so far given unphysical results, namely, the 
lattice starts to fill with baryons at a chemical potential well below the expected 
value of one-third the baryon mass. It seems plausible that this happens because 
the fi = ensemble does not overlap sufficiently with the finite-density states of 
interest, and so the true effects of quark loops will only be seen at exponentially 
large statistics @. 

In this paper we look at an alternative: evaluating the partition function 
at imaginary chemical potential, for which the measure remains positive, and 
standard Monte-Carlo methods apply. The canonical partition functions can then 
be obtained by a Fourier transform j|, Since the dominant source of errors is now 
the Fourier transform rather than poor overlap of the measure, it seems worthwhile 
to explore imaginary chemical potential as an alternative to the Glasgow method. 

An outline of the paper is as follows. We give criteria that a theory should 
satisfy in order for Monte-Carlo simulations at finite density to be feasible. We 
describe a toy model where even-odd effects effects become visible. We find some 
interesting examples (e.g. QCD at finite isospin density) where lattice simulations 
are possible. As a feasibility study we perform Monte-Carlo simulations for the 
two-dimensional Hubbard model with imaginary chemical potential, and find that 
it is indeed possible to obtain the canonical partition functions at low particle 
number. At the rather high temperature and low interaction strength that we 
study, we see no sign of electron pairing. 

2 Chemical potential and positivity of the measure 

Consider a generic system of fermions ip and bosons <ft, where the fermion Lagrange 
density is tfjM(<f))ifj. On integrating out the fermions, the partition function 
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becomes 

' V(j) e- 5 ^ w detM(</>). (2.1) 



In order to perform Monte-Carlo simulations it is necessary that the measure be 
nonnegative, so either we have to restrict ourselves to the cases where det M ^ 0, 
or to treat det M as an observable. The latter option is usually not viable, as 
det M tends to be a rapidly varying function of (f). We will discuss it again at the 
end of this section. 

To guarantee that the measure is positive, we must generally have an even 
number of flavors, for each of which detM is real (but not necessarily positive). 
One situation where det M is real is when there exists an invertible operator P 
such that 

M ] = PMP- 1 . (2.2) 

For a Wilson lattice fermion at zero chemical potential this relation holds, with 
P = 75, so any even number of flavors can be simulated by Monte-Carlo. With 
real chemical potential ( |2.2|) breaks down, but with imaginary chemical potential 
it is valid, and again simulations are possible for an even number of flavors. 

There are more exotic situations where fl2.2p holds. For example, consider two- 
flavor QCD with a finite density of isospin. In this case M has a block-diagonal 
structure 

MW=( L M ° Y (2.3) 



L{-fj) 

where \x is the chemical potential for the isospin, and L{p) is a Dirac operator 
for one flavor with chemical potential fi. L(fi) satisfies L(ji)' = 7s£(— /i)7s, hence 
Eq. Q2.2D is satisfied by setting 

n Y (2-4) 



. 2/5 . 

Here det M(/i) = | det L(fi)\ 2 ^ 0. More generally, consider QCD with Nf flavors. 
It has a global vector-like symmetry G = U(1)b x SU(Nf), where U(1)b is the 
baryon number. One may consider a nonzero chemical potential coupled to any 



generator T of G. Then Eq. (2.2) is satisfied for some choice of P if and only if 
the nonzero eigenvalues of T come in pairs A, —A. Thus det M is real for QCD 
at nonzero isospin density, but not for nonzero hypercharge density, or baryon 
number density. 

Another case where det M is real even in the presence of the chemical potential 
is when there exists an invertible operator Q such that 

M* = QMQ- 1 . (2.5) 
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Examples of this sort are afforded by models with four-fermion interactions, like 
the Hubbard model and the Gross-Neveu model. There M is a real operator, so 
det M is real too. Other examples are gauge theories with "quarks" in the real or 
pseudoreal representation of the gauge group. Thus det M is real for SU(2) with 
quarks in the fundamental representation, or for SU(3) with quarks in the adjoint, 
even when the chemical potential is nonzero. 

In some cases with det M real but not positive it appears that one can perform 
simulations by treating the sign of det M as an observable. The Hubbard model 
can be treated in this way far below half-filling |J. This is also the case for the 
Gross-Neveu model at nonzero chemical potential 0. Ref. jTj further argues that 
that det M is nonnegative for most of the configurations, so one can simply replace 
detM with |detM|. In QCD the phase of the determinant contains important 
physical information, but calculations have been performed without it [|J. 



3 Imaginary chemical potential 

The fact that the fermion determinant for QCD is real in the presence of an 
imaginary chemical potential \i — iv makes this an attractive option for exploring 
finite quark density. Simulations with an imaginary chemical potential are more or 
less equivalent to simulating a canonical ensemble. Indeed, the partition function 
for imaginary chemical potential 

Z(iv) = Tr e - pA e iPvft , (3.1) 

which is a periodic function of v with period 2ir//3, is the Fourier transform of the 
canonical partition function 



r-27r//3 

Z(N) = I dvZ{iv)e- iPvN . (3.2) 



o 



In principle, one can compute Z(iv) on a lattice as a function of v, and then use 
Eq. fl3.2p to obtain the canonical partition function. In practice, this method can 
work only for low enough N, because for large N the integrand in Eq. (|3.2| ) is a 
rapidly oscillating function, and the error of the numerical integration will grow 
exponentially with N. The method fails completely in the thermodynamic limit 
N — > oo. This need not discourage us, however, because in lattice simulations one 
is always working in a finite and rather small volume. The real question is how high 
we can push N before the numerical integration in Eq. ( |3.2|) becomes undoable. 
We will consider the two-dimensional Hubbard model as a testing ground for this 
approach. Related work has been performed in Ref. but without using the 
freedom to simulate at any v (see below). 
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Having a positive measure is not the end of the story. In practice we want to 
be able to use importance sampling to calculate Z(iv) with reasonable accuracy. 
To this end we rewrite Z{iv) in the following form: 



Now we treat the ratio of determinants as an observable, and the rest as the 
measure. A natural worry is that the ratio of determinants could be a rapidly 
varying function of the bosonic fields 0, which would make Monte-Carlo simulations 
unfeasible. This is what happens for real chemical potential. However, with an 
imaginary chemical potential this does not occur. The absolute value of the 
observable could still be a rapidly varying function. This has to be decided on 
a case-by-case basis. We will see that for the two-dimensional Hubbard model the 
ratio of determinants is a slowly varying function of <j), hence Z{iv) is computable. 

Using (|3.3| ) we can calculate the partition function for a range of v around a 
reference value vq. It is clear that we can cover the range u = . . . 2tt//3 with a set 
of "patches" each centered on a different uq. We can use as many patches as are 
required, so the measure will always overlap arbitrarily well with the observables. 

Some qualitative features of the system can be inferred from the knowledge 
of Z(iv) alone, without performing the Fourier transform. For example, consider 
a model of interacting fermions on a lattice (the example of the Hubbard model 
will be discussed further below). At some filling fraction the system may undergo 
a phase transition to a BCS superconducting phase. In that phase, the system 
will be populated with Cooper pairs, so the partition function will be dominated 
by sectors in which the charge is a multiple of 2. This will be clearly visible in 
Z{iv\ which will not only be perodic with period 2tt//3, but acquire a significant 
subharmonic at period 2tt/(2/3). This would be a signal that for some range of 
densities the energy of the system has a minimum when the number of electrons 
is a multiple of 2. 

This can be illustrated by a simple toy model containing a fermion with mass 
Mf and charge 1, and a boson with mass and charge 2. If we make less than 
Mf then, assuming some very weak charge-conserving interactions that establish 
equilibrium, states of even charge will be favored. 

The free energy is the sum of the fermion and boson contributions, 



where -Ff erm i on and -Fb OS on are the free energies of free fermions and bosons, 




(3.3) 



-ffermion -fboson 

(M 6 ,2//,T), 



(3.4) 
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Figure 1: Free energy as a function of imaginary chemical potential v (left) and 
log of the canonical partition function as a function of particle number (right) for 
our illustrative toy model ( |3.4|) , with j3M b = 1 and (3Mf = 5. 



respectively, 

J .fP- ± i n (i ± 2 cosh(/?/i)e- /3£ ^ + e~ 2 ^) : 
E{pf = p 2 + M 2 . 



fermion ' 
boson 



(3.5) 



The pairing of fermions into bosons is clearly visible in Z(N) (Fig. [T]), and 
also directly in Z(iu) which is not only perodic with period 2k / (3, but also 
approximately periodic with a smaller period 2ir/(2(3). This is a signal that the 
energy of a system has a minimum when the number of particles is a multiple of 
2. However, to infer the existence of an "unpairing" transition at high chemical 
potential, a visual inspection of the the plot of Z(iu) would not suffice: this 
information is encoded in the high-frequency behaviour of the Fourier transform 
of Z{iv). 



4 Hubbard model with imaginary chemical potential 

At densities away from half-filling, the single-flavor Hubbard model has a real 
but not necessarily positive fermionic determinant |J. Like QCD, the measure 
becomes positive at imaginary chemical potential. Since this model is of physical 
interest and also much less computationally demanding than QCD, it is interesting 
to use it to study the feasibility of performing simulations with imaginary chemical 
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potential. In fact, the model is so simple that for this initial investigation 
we were able to dispense with the usual hybrid Monte-Carlo algorithm |J for 
evaluating the fermion determinant, and perform the whole calculation with the 
computer mathematics tool "Mathematica" , using its "Det" function to calculate 
the fermionic determinants. 

Using the formulation described above (see ( |3.3|) ), we calculated the partition 
function as a function of imaginary chemical potential v. We used a 4 2 x 10 lattice 
with inverse temperature (3 = 1.5. (Lower statistics were also gathered on 4 2 x 20, 
to check that temporal discretization errors were under control.) The results are 
given in Fig. |^. By particle-hole symmetry, Ziiv) = Z(—iu), and Z has period 
2tt/P, so we only plot v = to tt//3. 

Three "patches" were used (see Sect. H), centered at (3v$ = 0, 7r/2, and ir. At 
the temperature we study, the error in Z(ii/)/Z(ivo) rises rapidly with \v — v$\, so 
it is crucial to use multiple patches to keep the statistical errors in Z{iv) under 
control. In turn, via the Fourier transform, this controls the errors in Z(N) for 
N > 0. In contrast, Ref. @] only used uq = 0, which is adequate only for the small 
volume, low temperature, and low particle number (N ^ 2) they studied. 

We then fitted Z{iv) to various ansatze, and Fourier transformed them to 
obtain the canonical partition function Z^ (Fig. [3]). Even using our very inefficient 
updating algorithm, we are easily able to to get accurate results up to N = 6. The 
error bars reflect the statistical errors in the Z{iv) data. We used fit functions of 
the form Z(iv) = exp(—au 2 ) x spline. 

It has been suggested that at some filling fraction the 2D Hubbard model may 
exhibit superconductivity, through pairing of the electrons to form Cooper pairs 
which then condense. Along the lines described in section |3], we would expect 
such pairing to manifest itself as an even-odd periodicity of Zn, leading to a 
characteristic extra bump in Z(iv). At the temperature and coupling that we 
studied, we see no such evidence of pairing. Obviously it would be very interesting 
to explore a wider range of parameters. This method might also be suitable for 
exploring the metal-insulator transition near half-filling, where the sign problem 
becomes particularly virulent ||. As noted in Ref. Q, however, it will become 
much harder to extract the Z^ at low temperatures, where the N = contribution 
dominates Z{iv). 

5 Conclusions 

Let us conclude with a few remarks concerning the possible utility of 
imaginary chemical potential in QCD. An imaginary chemical potential does not 
systematically bias the ensemble to large density, so that to pick out the effect of 
states of non-zero baryon number density one must rely on fluctuations (which, 
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Figure 2: Partition function for 2D Hubbard model as a function of imaginary 
chemical potential is, in units of 7r//3. This is on a 4 2 x 10 lattice, with (3 = 1.5, 
hopping term K — 1, and interaction strength [/ = 0.1, following the conventions 
of Creutz @. 



when they occur, are appropriately weighted by the action). These fluctuations 
will occur most readily when the temperature is high, and the gap in the baryon 
number channel is small. And even then one can only realistically hope, on the 
small lattices likely to be practical, to fluctuate to a few baryons. So a reasonable 
procedure would seem to be to start with a high temperature and work down, 
looking for qualitative chang function of temperature. In this way, one 

could realistically hope to use the methods described in this paper to study how 
properties of the quark-gluon plasma are affected by a net quark density. One could 
also study the deconfinement crossover, near which the baryons become light, and 
in particular locate the critical point in the T — /x plane predicted for two-flavor 
QCD PH| . All these phenomena are of immediate interest, since they will be 
explored in the next generation of heavy ion collision experiments. 

Another possibility is to work with large numbers of quark species, close to 
16, so that the theory is perturbative. Then there is no mass gap to baryons, so 
fluctuations are cheap, and also the contribution of interest, due to the quarks, is 
not swamped by gluons. In this case the cancellations may not be so bad even 
at real chemical potential, and the imaginary chemical potential approach should 
also work better, since the fluctuations of interest will occur frequently. 
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Figure 3: Log of canonical partition function Zjq for 2D Hubbard model, obtained 
by Fourier transform of Z{iv) (Fig. |^). 
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